##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--maf_file"), type = "character") ,
    make_option(c("--maf_msi_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    info_file <- paste(work_dir,"/config/STAD-useCombine.Sample.tsv",sep="")
    maf_file <- paste(work_dir,"/maf/All_ForMutBurden.extract.maf",sep="")
    maf_msi_file <- paste(work_dir,"/maf/All_ForMutBurden.extract.MSI.maf",sep="")
	images_path <- paste(work_dir,"/images/mutBurden",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_file <- opt$maf_file
maf_msi_file <- opt$maf_msi_file
info_file <- opt$info_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

col_im <- brewer.pal(9,"YlGnBu")[6:8]
names(col_im) <- c("IM(IGC)" , "IM(DGC)" , "IM(IGC_DGC)")

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################

dat_sample <- data.frame(fread( info_file ))
dat_maf <- data.frame(fread( maf_file ))
dat_maf_msi <- data.frame(fread( maf_msi_file ))

###########################################################################################
dat_maf$MS_Type <- "MSS"
dat_maf_msi$MS_Type <- "MSI"

dat_maf <- rbind(dat_maf , dat_maf_msi)
dat_maf <- subset( dat_maf , Variant_Classification %in% Variant_Types )

###########################################################################################
## 判断突变的共享情况
dat_sample <- data.frame(Patient = dat_sample$Patient , Class = dat_sample$Class , Tumor = dat_sample$Tumor , coverage_CDS = dat_sample$coverage_CDS)
dat_maf2 <- merge( dat_maf , dat_sample , by.x = "Tumor_Sample_Barcode" , by.y  = "Tumor" , all.x = T )
 
dat_maf2$Location <- paste( dat_maf2$Chromosome , dat_maf2$Start_position , 
    dat_maf2$Reference_Allele , dat_maf2$Tumor_Seq_Allele2 , sep = ":" )
dat_maf2$vid <- paste( dat_maf2$Location ,dat_maf2$Patient , sep = ":" )

pre <- "IM"
can <- c("IGC" , "DGC")

dat_maf2_share <- dat_maf2 %>%
group_by( vid ) %>%
summarize( Share = ifelse(length(which(Class %in%  pre))!=0 & length(which(Class %in% can ))!=0 , "Share" , ""  ))

dat_maf3 <- merge( dat_maf2 , dat_maf2_share , by = "vid" , all.x = "TRUE" )

###########################################################################################
## 计算WES的突变数量
exon_mut_share <- subset( dat_maf3 , Share == "Share"  ) %>%
group_by( Tumor_Sample_Barcode ) %>%
summarize( exonMutNum_Share = length(Start_position) )

exon_mut_private <- subset( dat_maf3 , Share == ""  ) %>%
group_by( Tumor_Sample_Barcode ) %>%
summarize( exonMutNum_Private = length(Start_position) )

exon_all <- merge( exon_mut_share , exon_mut_private , all.x = T , all.y = T )
exon_all[is.na(exon_all)] <- 0 

dat_plot <- merge( dat_sample , exon_all , by.x = "Tumor" , by.y = "Tumor_Sample_Barcode" , all.x = TRUE )

dat_plot$BurdenExon_share <- dat_plot$exonMutNum_Share/(dat_plot$coverage_CDS/1024/1024)
dat_plot$BurdenExon_private <- dat_plot$exonMutNum_Private/(dat_plot$coverage_CDS/1024/1024)

dat_plot$BurdenExon_share[is.na(dat_plot$BurdenExon_share)] <- 0
dat_plot$BurdenExon_private[is.na(dat_plot$BurdenExon_private)] <- 0

out_name <- paste0( images_path , "/mutBurden.share_private.tsv" )
write.table( dat_plot , out_name , row.names = F , sep = '\t' , quote = F  )